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AN EFFICIENT ALGORITHM USING MATRIX METHODS TO SOLVE 
WIND-TUNNEL FORCE -BALANCE EQUATIONS* 

By David L. Smith 
Langley Research Center 

SUMMARY 

An iterative procedure applying matrix methods to accomplish an efficient algorithm 
for automatic computer reduction of wind-tunnel force -balance data has been developed. 
Balance equations are expressed in a matrix form that is convenient for storing balance 
sensitivities and interaction coefficient values for online or offline batch data reduction. 
The convergence of the iterative values to a unique solution of this system of equations 
is investigated, and it is shown that for balances which satisfy the criteria discussed, 
this type of solution does occur. Methods for making sensitivity adjustments and initial 
load effect considerations in wind-tunnel applications are also discussed, and the logic 
for determining the convergence accuracy limits for the iterative solution is given. 

This more efficient data reduction program is compared with the technique pres- 
ently in use at the NASA Langley Research Center, and computational times on the order 
of one-third or less are demonstrated by use of this new program. 

INTRODUCTION 

Since aircraft and space vehicle motions depend on the forces and moments about 
the three orthogonal body axes, an extensive amount of wind-tunnel testing is devoted to 
measuring these quantities for given model configurations to enable the estimation of 
aerodynamic loads on full-scale vehicles in flight. The most commonly used method for 
measuring these forces and moments is by installing an internal strain-gage balance 
within a wind-tunnel model as illustrated in figure 1. The model is attached to this bal- 
ance and the forces and moments about the axes shown in figure 2 are transduced into 
electrical signals suitable for analog -to -digital conversion and subsequent data reduction 
or online evaluation where such equipment exists. 

*The information presented herein is largely based on a thesis entitled "The Appli- 
cation of Matrix Methods to Solving Wind-Tunnel Force-Balance Equations” submitted by 
the author to the Faculty of the School of Engineering and Applied Science of George 
Washington University in partial satisfaction of the requirements for the degree of 
Master of Science, December 1971. 



Wind-tunnel 

model 




Figure 2.- Force and moment axes with positive directions shown. 



The purpose of this investigation is to apply matrix methods to the force-balance 
equations in order to develop an efficient data reduction program which offers signifi- 
cantly fewer arithmetic operations and smaller computational times per data point. This 
program uses an iterative procedure to account for balance interactions and considers 
required sensitivity adjustments and initial load effects. A description of the balance 
data reduction is given and a technique is presented for determining from calibration 
data whether the iterative procedure will converge. The technique presently in use at 
the Langley Research Center is presented in appendix A. 


SYMBOLS 

Measurements are given in both SI and U.S. Customary Units. The measurements 
and calculations were made in U.S. Customary Units. The force and moment axes 
usually coincide with the body axes of a wind-tunnel model as shown in figure 2. 

c element of matrix M 

F^ axial force 

normal force 

Fy side force 

f generalized force or moment function 

J upper bound for Lipschitz’ constant 

K normalized coefficient 

k calibration coefficient 

M-^ rolling moment 

My pitching moment 

M^ yawing moment 

X generalized force or moment component 

a angle of attack 
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e nonlinear interaction correction 

8 meter indication 

k sensitivity 

Matrices: 

B positive full-scale balance design loads matrix (6 x 1) 

Cj first-order coefficient matrix (6 x 6) 

Cg nonlinear interaction coefficient matrix (6 x 21) 

E second -order interaction column matrix (6 x 1) 

F force and moment column matrix (6 x 1) 

X sensitivity diagonal matrix (6 x 6) 

M matrix product of (6 x 21) 

P force and moment product matrix (21x1) 

0 output column matrix (6 x 1) 

Subscripts: 

1 force or moment component considered 

j interacting load (see table I) 

k data point specified 

max maximum value 

min minimum value 

n number of iteration 
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o 


initial load effect 


u uncorrected value 


DEVELOPMENT OF MATRIX RELATIONS 


Background for Analysis 

Ideally, the output for each force or moment component measured by a balance 
should be affected only by a loading on that particular component. Experience shows, 
however, that a given component is often affected by loading another component. This 
effect is called an "interaction.” Interactions are classified as either linear or non- 
linear, depending on whether they are related to a single component's load or to exponen- 
tial powers and combinations of the components being loaded. Linear interactions result 
from machining tolerances, strain-gage positioning tolerances, variations in strain-gage 
properties, electrical circuitry, or Poisson's effect. Nonlinear interactions are attrib- 
uted to deflections of the strain-gage beams (ref. 1). 

Consider the general case of a balance designed to measure three perpendicular 
forces and three moments about the axes of these forces. The output of each component 
is a function of all six forces and moments due to presence of interactions and can be 
expressed as a polynomial of the form (ref. 2): 


Ot = k. * F, t + k. 0 F a + k. 0 M Tr + . . . + k, „F„ + k. n F* + k. Q F XT F 


i - *i,r N T 1,2 A A T i,3 Y 


1,6^ T i,7*N i,8 N A 


+ k i,9 F N M Y + • • • + k i,27 F Y 2 + k i,28 F N 3 + k i,29 F N 2r A + * ‘ * W 


The linear interaction coefficients for this case would be k. j, k^ • • •> g or the 
coefficients of the first-order terms of equation (1), except for the k^j which is the 
inverse of the ith component's sensitivity. Nonlinear interaction coefficients would be 
k. 7 , k. ft , . . ., or coefficients of second degree and higher order terms. In practice, 
third and higher order interaction terms are negligible, and second-order terms are 
generally small compared to the linear terms. (See refs. 1 and 2.) 

To facilitate force -balance data reduction, equation (1) is divided through by k* i 

i,i 

or "normalized,” with third and higher order terms neglected, yielding (ref. 2): 


( X i)u = K M F N + K i,2 F A + K i,3 M Y + • * * + K i,6 F Y 

+ K i,7 F N 2 + K i,8 F N F A + K i,9 F N M Y + • • • + K i,27 F Y 


(2a) 
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where 


k. - 


K- • = — ^ = Normalized interaction coefficients when i / j 

*>J ki j 




with 


K; i = = ith component sensitivity when i = j 

1,1 k i,i 


and 


(Xj) = k- :6- = Uncorrected force or moment on ith component 
' 'u U 1 1 

For a typical balance load, for example, normal force, (Xj) u = and 

kf i 

K, , = — *— = 1 which results in the following form of equation (2): 

1,1 k l,l 


F N " ( F n) u " ( K 1,2 F A + K 1,3 M Y + • • • + K 1,6 F Y 

+ K 1,7 F N 2 + K 1,8 F N F A + K 1,9 F N M Y + ' * ’ + K 1,27 F Y 2 ) 


(2b) 


In this form the interaction coefficients are expressed in terms of the apparent load on 
the ith component per unit of the jth loading. 


Assumptions 

In order to solve the system of six force-balance equations represented by equa- 
tions (2) for the aerodynamic loads acting on a wind-tunnel model, the interaction coef- 
ficients and sensitivity constants used in acquiring the data must be known. It will be 
assumed that these constants are available from the calibration of the balance, and that 
the sensitivities have been adjusted to the actual values in the tunnel installation. It will 
also be assumed that there are no initial load effects to account for at this time. Methods 
for including both the balance sensitivity adjustments and the initial load effects will be 
considered under a subsequent heading. 


Matrix Relations 

By defining the following matrices, equation (2a) can be expressed as a matrix 
relation where the subscripts denote the loads indicated in table I: 
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TABLE I.- LOAD CORRESPONDING TO A GIVEN SUBSCRIPT 


Load 

Load denoted by 
subscript - 

Load 

Load denoted by 
subscript - 

f n 

i 

f a m x 

15 

f a 

2 

f a m z 

16 

M y 

3 

f a f y 

17 

M x 

4 

M y 2 

18 

M z 

5 

M y M x 

19 

f y 

6 

M y M z 

20 

F 2 
r N 

7 

M y F y 

21 

f n f a 

8 

M 2 
M X 

22 

f n m y 

9 

M X M Z 

23 

f n m x 

10 

m x f y 

24 

f n m z 

11 

M 2 

z 

25 

f n f y 

12 

m z f y 

26 

f a 2 

13 

F 2 
r Y 

27 

f a m y 

14 




The output column matrix: 



(i = 1, 2, . . 6) 
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The force and moment column matrix: 


F = 




(i = 1, 2, . 


The second-order force and moment product column matrix: 


P = 


F 2 
F N F A 

f n m y 



/i = 1, 2, . . 6 

\j = i> i + 1, • • •* 



The 6x6 sensitivity diagonal matrix: 


“M 

0 

0 

0 

0 

0 


0 

* 2,2 

0 

0 

0 

0 


0 

0 

* 3,3 

0 

0 

0 


0 

0 

0 

* 4,4 

0 

0 


0 

0 

0 

0 

* 5,5 

0 


0 

0 

0 

0 

0 

* 6,6 
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The 6x6 first-order interaction coefficient matrix: 



’ 1 K l,2 

K l,3 

k m 

K l,5 

K l,6~ 


*2,1 1 

“2,3 

K 2,4 

“2,5 

K 2,6 

0 

11 

K 3,l • 

* 

■ 

* 

• 


K 6,l K 6,2 

K 6,3 

K 6,4 

K 6,5 

1 


The 6x21 second-order interaction coefficient matrix: 


K l,8 K l,9 
K 2,8 K 2,9 

K 6,8 K 6,9 

By use of these matrices, the six equations represented by equation (2) can be expressed 
as follows: 

itO^CjF+CgP (3) 

By use of matrix algebra (ref. 3), equation (3) is readily solved for F by subtracting 
C 2 P from each side and premultiplying each term by the matrix inverse of to 

obtain 

F = - Cj^CgP (4) 

In this form, equation (4) is very convenient for using iterative procedures to solve 
for F, or the forces and moments acting on the model. An iterative procedure is the 
most practical method of solving this equation because the second-order interaction cor- 
rections Cj”*C 2 P are functions of the elements of F. 

Since the previously defined coefficient matrices are made up of constants deter- 
mined from calibration data, and the product of C^”lC 2 can be calculated and 

stored in this form for subsequent data reduction. Also, if a balance is designed to mea- 
sure less than six components, the coefficient matrices can be accordingly reduced in 
size before these calculations are made. Carrying out these steps prior to actually 
reducing tunnel data greatly increases the efficiency of the data reduction program. 
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Iterative Procedures 


Upon examination, equation (4) tends to appear cumbersome or awkward to solve 
iteratively as each term on the right-hand side is a product of three matrices. Further 
examination shows this is not the case; actually, it is in a rather convenient form for the 
data reduction program. The product £0, or the engineering unit conversion, is cal- 
culated prior to the iteration stage of the data reduction and is called the "uncorrected" 
load, F u . Also, the product of Cj'iCg is calculated from the calibration data and is 
stored as a 6 x 21 matrix, M. Consequently, equation (6) can be expressed as follows: 

F=C 1 ’ 1 F U -MP (5) 

For a given data point G^, Fj = C 1 " 1 F U is directly calculable and is a very good 
approximation of F since Fj contains all first-order interaction corrections and 
because of the relatively small effects of second-order balance interactions. For this 
reason, the elements of Fj are used as the values of the forces and moments necessary 
for calculating the elements of the first approximation of the second-order matrix. Note 
that Fj^ = C^KjO is dependent only on calibration constants and on the meter indica- 
tions for the particular data point being reduced. These linear terms are directly cal- 
culated and require no iterating for their evaluation. Only the second-order interaction 
terms are iterated until the procedure converges. 

Iteration of the second-order terms is accomplished as follows: The first approx- 
imation of the second-order matrix Pj is premultiplied by M and generates the 
second-order interaction correction matrix The absolute value of each E 1 ele- 

ment | e A | is then compared with the required accuracy of convergence for each balance 
component. If each j e ^ | is less than the given convergence limit, which can be speci- 
fied or is calculated based on the balance sensitivities, the data reduction is complete 
for that data point with F = Fj - MPj. However, if one or more of the |e i | is greater 
than these convergence accuracies, equation (5) is reiterated as follows: 


Fn-V Fu-Vi 


(6) 


The quantity F n is then used to reevaluate P n which is again premultiplied by M to 
determine more nearly exact values of the second-order balance interactions E n . The 
column matrix E n is then compared with E ^ and if corresponding elements agree 
to within the convergence limits, the force and moment matrix may be expressed as 


F = Ci 


-1 


F - E 
r u n 


(V) 


Equation (6) is reiterated until convergence occurs. 


10 



Convergence 


In any iterative solution such as the one described, certain questions must be con- 
sidered such as whether F n will always converge, and whether its limit is a unique 
solution of the given equations. Henrici (ref. 4) considers these questions for the general 
case of iterating systems of nonlinear equations. A criterion for proving that equation (5) 
will converge is given and discussed in appendix B. Theorems given therein not only 
prove that this iterative procedure will converge but also show that successive iterations 
approach a unique solution in the region of the design loads for balances which satisfy the 
given conditions and have an upper bound for Lipschitz T constant J of less than 1 
where 



A computer subroutine program which calculates the value of J is given in appendix B. 

CONSIDERATIONS FOR REDUCING BALANCE DATA 
Sensitivity Adjustments 

The reduction of force -balance data requires each force and moment component 
calibration sensitivity to be adjusted to the actual values for the wind-tunnel installation. 
To accomplish this sensitivity adjustment, the same apparent loading is applied to each 
component at the calibration facility and at the wind-tunnel installation. The ratio of the 
output at calibration to the output at the tunnel installation for this common apparent 
load is then multiplied by the corresponding calibration sensitivity and yields the tunnel 
sensitivity as follows: 



tunnel 



calibration 


calibration ^ 
^tunnel ) 


( 9 ) 


An efficient method for making this sensitivity adjustment is to store the calibration sen- 
sitivities and apparent load outputs with the balance interaction coefficient values. Then, 
adjusted tunnel sensitivities can readily be computed and assigned to the proper locations 
in the X matrix by supplying balance output values for the same apparent load in the 
tunnel installation, and performing the indicated ratios in the data reduction program. 
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Initial Loads 


Another important consideration in force-balance data reduction stems from the 
fact that the balance equations (2) are nonlinear. For this reason, tunnel data must be 
related to the same origin as calibration data, or zero output for zero loads on all com- 
ponents, as shown by the solid line in figure 3. Typically, meter readings at wind-off 
zero-angle-of-attack conditions are taken as the zero load values or as the origin of the 
data. However, initial loads such as model weight cause the balance output to be located 
off the calibration origin, for example, at point A in figure 3 where the prime indicates 
the tunnel axis system. Ignoring initial load effects is essentially the same as assuming 
that the balance is performing according to the dashed curve or the calibration curve 
shifted to the new origin. The balance is actually performing according to the solid cal- 
ibration curve. Therefore, the data origin must be shifted to correspond with the cali- 
bration origin for each data point. 
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A convenient method for reducing balance data with initial loads is to translate the 
axes to the system used for calibration, eliminate balance interaction effects, and then 
translate the axes back to the set used in acquiring the data. Note in figure 3, that 0 O 
and X 0 are determined prior to taking aerodynamic data, and 0' is recorded for each 
data point during a wind-tunnel test. This observation suggests the following axes 
translation: 


0= e' + e 0 (io) 

Substituting this value into the balance equation allows the calculation of X, from which 
X' is readily determined by the following translation back to the primed axes: 

X' = X-X 0 (11) 

This method is readily extended for a six-component balance as shown in the following 
matrix relations: 

0 = 0' + 0 O (12) 

This output matrix is then substituted into equation (5) and is solved iteratively as 
described for F, from which 

F* = F - F 0 (13) 

This method of translating axes to include initial load corrections in balance data reduc- 
tion has been used with the iterative procedure discussed previously. The only arithmetic 
operations required for these axes transfers are six additions before iterating the balance 
equations, and six subtractions following the iterations. 

An alternate method for considering initial loads by reevaluating balance interaction 
coefficients to account for these axes translations is described in reference 5. 

Convergence Limits 

For any iterative solution, an accuracy or convergence limit must be specified. 

This limit can be an absolute value, as presently used with balance data at Langley, a 
percentage of the solution itself, a percentage of the maximum range of the solution, or a 
percentage of the resolution of the data acquisition system. Because of the limitation of 
the recording system resolution, the minimum detectable increment of each component is 
equal to its sensitivity times 1 count, or 
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( AX i)min = *i,i X 


1 = K- 


M 


(14) 


Consequently, the convergence limit of one -tenth the value of this minimum detectable 
increment, or k ± i/10, for each force and moment component is used in this data reduc- 
tion program. This criterion will cause to be negligible any systematic errors that may 
result because of the convergence accuracy. 


Data Reduction Program 

The iterative procedure and other related topics that must be considered for balance 
data reduction have been utilized in developing a FORTRAN program for the Control Data 
6600 digital computer complex at the Langley Research Center. This program is listed 
and described in appendix C. 


COMPARISON OF PRESENT DATA REDUCTION METHOD 
WITH MATRIX METHODS 


Logic 

The logic of the matrix method developed in this paper and of the technique pres- 
ently in use at the Langley Research Center (appendix A) is very similar in many ways. 
Both methods apply an iterative solution of the form 

F n=l(f„-l) (15) 

where _f denotes a column vector. (See appendix B.) The presently used method iter- 
ates each force and moment component individually and updates or recalculates the 
second-order products between each components iteration. The matrix method, however, 
iterates the column matrix F n or updates all force and moment components before the 
second-order product matrix is recalculated. Also, the first approximation of these two 
methods is determined differently. The present method uses the products of the sensitiv- 
ities times the balance outputs or k* -0. and iterates both first- and second-order 
interaction terms, but the matrix method uses C as the first approximation and 
consequently must iterate only the second-order terms. 

Computation Time Requirements 

Because of the differences described, the number of arithmetic operations required 
per data point by using the matrix method is considerably reduced and, as a result, corre- 
sponding decreases in the computation time for each data point are obtained. Table II 
gives a comparison of the deviation of the first approximation iterated in the balance 
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equations for the two methods. These values are significantly closer to the calculated 
solutions for each balance component in the matrix method. Consequently, fewer itera- 
tions are required for the data points in using the matrix method as shown in three of the 
four cases presented in tables III and IV. There are also significant differences in the 
number of arithmetic operations for the two methods due not only to the few iterations 
but also to the facts that load combinations are updated after iterating all components and 
only second-order terms must be iterated by using the matrix method. This reduction 
in the number of arithmetic operations results in computation times on the order of one- 
third or less for the matrix method over the present technique. Table V shows consid- 
erable reductions in the number of arithmetic operations and computation times even 
when both methods are driven through the same number of iterations. These observations 
are especially noticeable for balances with less than six components because of the 
"collapsing" of the coefficient and data matrices as discussed previously. The matrix 
method is thereby significantly more efficient than the present technique in which the 
coefficients for components not measured are set equal to zero and the arithmetic oper- 
ations are performed for all components and, as a result, there is the same number of 
calculations per iteration for any number of measured force and moment components. 


TABLE II.- DEVIATION OF FIRST APPROXIMATION FROM SOLUTION 
FOR BOTH METHODS OF DATA REDUCTION 


Balance 

component 

Converged 

iterative 

solution 

First approximation 
values 

Deviation of 
approximation from 
iterative solution 

Present 

method 

Matrix 

method 

Present 

method 

Matrix 

method 

F N , N (lb) 

350.5 (78.8) 

343.0 (77.1) 

350.1 (78.7) 

7.5 (1.7) 

0.4 (0.1) 

F A , N (lb) 

53.8 (12.1) 

67.6 (15.2) 

55.6 (12.5) 

13.8 (3.1) 

1.8 (0.4) 

My, N-m (in-lb) . . 

10.1 (89.6) 

9.3 (82.0) 

10.0 (88.9) 

0.8 (7.6) 

0.1 (0.7) 

M^, N-m (in-lb) . . 

3.03 (26.8) 

3.17 (28.1) 

3.03 (26.8) 

0.14 (1.3) 

0 (0) 

M z , N-m (in-lb) . . 

5.65 (50.0) 

5.62 (49.7) 

5.74 (50.8) 

0.03 (0.3) 

0.09 (0.8) 

F y , N (lb) 

130 (29.2) 

166 (37.3) 

131 (29.4) 

36 (8.1) 

1 (0.2) 
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TABLE m.- COMPARISON OF MATRIX METHOD WITH PRESENT TECHNIQUE 
FOR TYPICAL DATA WITH NO INITIAL LOAD TRANSLATIONS 


Number of 
balance 
components 
reduced 

Number of 
iterations 

Approximate 
number of 
arithmetic 
operations 

Time required for 
iterations, msec 

Present 

method 

Matrix 

method 

Present 

method 

Matrix 

method 

Present 

method 

Matrix 

method 

3 

6 

4 

1728 

105 

16 

2 

4 

3 

3 

864 

166 

12 

2 

5 

3 

2 

864 

205 

10 

4 

6 

5 

2 

1440 


18 

6 


TABLE IV.- COMPARISON OF MATRIX METHOD WITH PRESENT 
TECHNIQUE FOR TYPICAL DATA WITH INITIAL 
LOAD TRANSLATIONS REQUIRED 


Number of 
balance 
components 
reduced 

Number of 
iterations 

Approximate 
number of 
arithmetic 
operations 

Time required for 
iterations, msec 

Present 

method 

Matrix 

method 

Present 

method 

Matrix 

method 

Present 

method 

Matrix 

method 

3 

6 

4 

1740 

111 

16 

4 

4 

3 

3 

876 

174 

12 

4 

5 

3 

2 

876 

215 

12 

4 

6 

5 

2 

1452 

342 

18 

4 
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TABLE V.- COMPARISON OF MATRIX METHOD WITH PRESENT 


TECHNIQUE FOR TYPICAL DATA POINTS WITH THE SAME 
NUMBER OF ITERATIONS FOR BOTH METHODS 


Number of 
balance 
components 
reduced 

Number of 
iterations 

Approximate 
number of 
arithmetic 
operations 

Time required for 
iterations, msec 

Present 

method 

Matrix 

method 

Present 

method 

Matrix 

method 

Present 

method 

Matrix 

method 

3 

7 

7 

2016 

171 

18 

4 

4 

4 

4 

1152 

208 

12 

4 

5 

3 

3 

864 

295 

10 

6 

6 

4 

4 

1152 

612 

14 
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CONCLUDING REMARKS 

The reduction of wind-tunnel force-balance data by applying matrix methods to an 
iterative solution of the balance equations has been presented, and it has been demon- 
strated that this method is a significant improvement over the presently used method at 
Langley Research Center. The convergence of this iterative solution was considered, 
and it was shown that for balance equations which satisfy the conditions specified, this 
method would converge to a unique solution within the range of the design loads of the 
balance. A technique was also presented to determine whether the balance equations 
satisfy these conditions based on the calibration data for the balance. Considerations of 
sensitivity adjustments and initial load effects were discussed and methods for making 
these corrections were given. 

This matrix method has been developed with the assumption that the third and 
higher order balance interactions are negligible. If the case arises in which such inter- 
actions must be considered, this calculation can be readily accomplished with these 
methods by adding the load combination (s) which produce the third or higher order inter- 
action to the force and moment product column matrix and including the appropriate 
coefficients on each row of the nonlinear interaction coefficient matrix. These changes 
would, of course, necessitate changing the dimensions of these arrays in the computer 
programs given herein. 
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The efficiency of this force-balance data reduction algorithm resulting from 
applying matrix methods makes it particularly useful for real-time display and control 
calculations by smaller online computers as well as beneficial for offline batch data 
reduction subsequent to the wind-tunnel test runs. Computational times of one-third or 
less than those required for the presently used technique are demonstrated by this matrix 
methods algorithm. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., June 7, 1972. 
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APPENDIX A 


A PRESENT METHOD OF FORCE -BALANCE DATA REDUCTION 

The present method used to reduce force-balance data at the Langley Research 
Center involves an iterative solution of the six balance equations represented as follows: 

( X i)u = K i, 1 F N + K i,2 F A + K i,3 M Y + ■ • • + K i ? 6 F Y 

+ K i,7 F N 2 + K i,8 F N F A + K i,9 F N M Y + • • • + K i,27 F Y ( A1 ) 

To solve these equations, the calibration sensitivity and all first- and second-order inter- 
action coefficients must be known for each force and moment component. 

After each force and moment component calibration sensitivity is adjusted to its 
actual value for the wind-tunnel installation, data reduction is accomplished through use 
of a computer subroutine program based on iterating equation (Al). The first approxi- 
mation of each force and moment is obtained from 



K- 


1,1 1 


(A2) 


where the prime indicates that the tunnel sensitivity adjustment has been made. These 
first approximations of the aerodynamic loads are added to the initial load values and 
then substituted into equation (Al) to calculate a first approximate value of the interaction 
correction for each force and moment component. These correction values and the initial 
loads are then subtracted from the approximations of equation (A2), and the results become 
the second approximations of the aerodynamic loads. These second approximations are 
then added to the initial loads and reiterated into equation (Al), from which a second 
approximation to the interaction corrections is determined. The first and second inter- 
action correction approximations are compared for each balance component and if they 
agree with a specified tolerance, then the latter corrections are subtracted from the 
force and moment approximations and the balance data reduction is complete. If these 
two approximations do not agree within the given tolerance, then the latter interaction 
corrections are subtracted from the force and moment approximations from equation (A2) 
and these values are reiterated into equation (Al) until convergence occurs for all balance 
components. 
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APPENDIX B 


CONVERGENCE OF ITERATIVE SOLUTIONS 

In order to develop a criterion for the convergence of the system of nonlinear bal- 
ance equations, it will be convenient to use vector notation. The coordinates of the point 
(F n , F^, m y’ M X’ M Z’ F y) can be re P resented by the column vector F or 


F = 


f n 


f l( F N’ F A’ M Y’ M X’ M Z’ F y) 

f a 

— 

^2 ( F N * F A’ M Y’ M X’ M Z’ F y) 

f y_ 


f 6( F N’ F A’ M Y’ M X’ M Z’ F y)_ 


(Bl) 


It is also convenient to denote a column vector with elements of f j, f 2 * 
Equation (Bl) can thusly be written as follows: 

F = f(F) 


f 6 as 1(F). 
(B2) 


By using this notation, the following theorem given by Henrici (ref. 4) can be applied to 
the force-balance equations or to equation (5): 

Theorem. Let R denote the region with limits aj and b^ 

a l= F N= b l 
a 2= F A= b 2 
a 3 = M Y = b 3 
a 4 = M X = b 4 
a 5= M Z= b 5 


a 6 = F Y = b 6 
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APPENDIX B - Continued 


and let the functions _f satisfy the following conditions: 

(a) f j, fg, . . fg are defined and continuous on R. 

(b) For each FeR, the point fj(F), f 2 (F), . . f g (F) also lies in R. 

(c) There exists a constant L < 1 such that for any two points Fj and F 2 
in R, the following inequality holds: 


|i( F l)-i( F 2)|| SL 



(B3) 


where the double bars denote the Euclidean norm. Then the following statements are 
true: 


(a) Equation (B2) has precisely one solution S in R. 

(b) For any choice of F Q in R, the limit of the iterative solution described or 
F n = l(F n _i) is defined and converges to the unique solution S. 

(c) For any n = 1, 2, . . ., the following inequality holds: 



F l- 



(B4) 


It can easily be shown that the expressions of f^, f 2 , . . fg satisfy conditions (a) 
and (b). Let the region R be bounded by 1.5 times the minimum and maximum loads 
for which a balance is designed to measure. The 1.5 factor is necessary as interaction 
effects can cause the first approximations to be outside of the design load region. Now 
consider f^ expanded as follows: 


'i( f n- f a- m y- m x> m z- f yH f n), + c i, 7 ( f n) 2 


c 1,8 F N F A 


+ c 


1,27 ( F y) 


(B5) 


where 


( F n)j 1 


is the first approximation of normal force or the first element of 


is con- 


C 1 " 1 F U . For a given data point 0^, (F n ) is constant. It is obvious that fj i 

tinuous in R as are f 2 , . . ., fg (ref. 6), and by virtue of equations (Bl) and (B5), con- 
dition (b) also is satisfied. 


In order to establish that condition (c) is satisfied by the second-degree expressions 
f j, f 2 , . . ., fg, the Lipschitz constant L must be evaluated or a maximum for its value 
must be established. Henrici (ref. 4) has developed a criterion for determining the bound 
of the Lipschitz constant, which is given in the following theorem: 
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APPENDIX B - Continued 


Let the functions fj, f 2 , . . fg have continuous partial derivatives in the 
region R as defined. Then, the inequality (B3) holds with L = J, where 



/ i = 1, 2 6\ (B6) 

U = 1, 2, . . 6 / 


The limiting value of J is calculated rather straightforwardly by taking the par- 
tial derivatives of the balance equations and evaluating the maximum possible value of 

9 f l 

each term as follows for — — : 

aF N 


af 

8F 


— - 2Cj 7 F n + c 18 F a + Cj 9 M y + Cj 10 M X + C 1 11 M Z + c l, 12 F Y 
N 


(B7) 


Each term on the right-hand side of equation (B7) is evaluated at 1.5 times the maximum 
design load, and the absolute values of these products are summed. This method of 
evaluation is carried out for each partial derivative in equation (B6) and results in an 
upper bound for J or 



A computer subroutine program which calculates the upper bound for J is listed. 
This program assumes that the C 1 “ 1 C 2 product array is stored in M and that the 
design loads are stored in a one-dimension array B. The maximums of the partial 
derivatives are computed as indicated above and stored in the 6 X 6 array A from 
which the upper bound for J is calculated in accordance with equation (B8). The 
1.5 factor is not applied in this subroutine but should be considered when interpreting the 
result of this evaluation. In practice this factor can be varied depending on the size of 
the interactions on the balance. 

Because this upper limit of the Lipschitz constant is dependent only on the inter- 
action coefficients and the design loads of a balance, it can be evaluated prior to the use 
of a balance in a wind-tunnel application. It is convenient to determine the bound of this 
constant at the same time the interaction coefficients are evaluated. 


22 



APPENDIX B - Concluded 


c 

c 

c 

c 

c 


10 


SUBROUTINE LPSHZ(M,B, IBAL* IDATE) 

DIMENSION M(6*21> B<6) 

THE VALUES OF REQUIRED PARTIAL DERIVATIVES WILL BE EVALUATED FROM 
THE DESIGN LOADS STORED IN B AND THE SECOND-ORDER INTERACTIONS 
STORED IN M. THE PARjIAl DERIVATIVES WILL BE STORED IN A AND THE 
LIPSCHITZ CONSTANT STORED IN AL I P * IBAL AND IDATE ARE THE 
BALANCE NAME AND CALIBRATION DATE RESPECTIVELY IN DISPLAY CODE* 

DO 10 I “1 * 6 

( ABS (2**M < I • j )*B ( 1 ) ) +AB S ( M (I*2)*B(2>)+ABS(M(I*3)*B<3))+ 

ABS ( M ( I *4 )*B (4 ) > +AB S ( M ( 1 *5 )*B (5 ) )+ABS ( M ( I * 6 ) *B < 6 ) ) ) 
(ABF ( 2« *M ( I • 7 ) *B ( 2 ) ) + ABS ( M ( I *2 )*9 ( 1 > )+ABS <M(I*8>*B<3>)+ 

M ( I *9 )#B (4 > )+ABS (M ( I • 1 0 >*B < 1 ) ) +ABS (M ( I * 1 1 )BS(6) ) ) 

“ * 2 )*B( 3 ) )+ABS (M ( I • 3 )*B ( 1 ) > + ABS (M ( 1 * 8 )*B (2 ) > + 

ABS (M < I * 1 3 ) *B ( 4 ) > + ABS (M ( I * 1 4 > *B (5 ) ) + 
ABS ( M ( I * 15 )*B (6) > ) 

6 )*B (4 ) l + ABS (M { I * 4 )*B < 1 ) >+ABS (M ( I « 9 )*B (2 ) ) + 
ABS (M ( I * 1 3 >*B<3 > >+ABS(M ( I * 1 7 > *B ( 5 ) > + 
ABS ( M ( I * 18 ) *B (6 ) ) ) 

9 >*0 (5 ) ) + ABS (M { I « 5 )*B ( 1 ) ) + ABS ( M ( ! , 10 )*B (2) ) + 
ABS (M ( I • 1 4 )*B(3 ) ) + ABS ( M ( I « 1 7 ) *B ( 4 > ) + 
ABS(M( I *20 ) *8 ( 6 > ) ) 

( ABS (2**M ( I tol )*B (6 ) ) + ABS (M ( I *6 )*S9l ) ) +ABS ( BM 1*11) *SB2 ) ) + 

ABS (M ( I • 1 5 )*B ( 3 ) ) + ABS ( M ( I * 1 8 > *8 (4 ) > + 
ABS ( M ( I *20)*B<5) ) > 

CONTINUE 

CALCuALTlON OF LIPSCHITZ constant* ALIP 
AL I P»0 • 0 
DO 20 L*1 *36 
AL IP*AL I P+A (L )**2 


A < I * 1 ) * 

( ABS ( 2 * *M ( 


ABS ( M ( 

A < I * ? ) = 

( ABF ( 2 • *M ( 


ABS ( M< 

A ( I « ? ) * 

(ABF ( 2 **M ( 

A ( I ♦ 4 ) * 

( ABS (2**M ( 

A ( I « B ) * 

( ABS ( 2 • *M ( 

A < I *6 > » 

( ABS ( 2 • *M ( 


20 CONTINUE 

ALlP*SQRT < AL I P ) 

PRINT 15* IBAL* IDATE 

15 FORMAT ( 1 HI « ///////* BALANCE **A10* # DATE **A10) 

PRINT 16 

16 FORMAT(1HO** THE FOLLOWING ARRAy CONTAINS THE PARTIAL DERIVATIVES 

1 OF X(U) WITH RESPECT TO X( !>•*/* WHERE (I) DESIGNATES THE ROW AND 

2 <J) DESIGNATES THE COLUMN**///) 

PRINT 1 7 • A 

17 FORMAT ( 1HO*6(2X*F10«6>/) 

PRINT 18* ALIP 

IS FORMAT ( 1 HO* * THE LIPSCHITZ CONSTANT FOR THIS BALANCE IS LESS THAN 

S*«F1 0*5 ) 

RETURN 

END 
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APPENDIX C 


BALANCE DATA REDUCTION SUBROUTINE PROGRAMS 

The subroutine programs given in this appendix correct force -balance data for 
interaction effects by applying the matrix methods discussed. Subroutine CTRNL calcu- 
lates the initial load corrections necessary for determining second-order interactions on 
a multicomponent balance. Subroutine CINTR then corrects balance data for both first - 
and second-order interactions, considering initial load effects where required. Provi- 
sions are also included to account for one discontinuous interaction term, that is, an 
interaction coefficient for which the value depends on whether a particular component's 
load is positive or negative. 

A flow chart of the subroutine CINTR follows. The listings of the two subroutine 
programs CTRNL and CINTR along with the required matrix operations subroutines are 
given with pertinent comments after the flow chart. 
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c 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE C I N TR ( FU ,F l , E l , L I ST , F , I ER ) 


********** 4 4* ************************* ***************************** 


* * 

* SUBROUTINE * 

* CINTR * 

* * 

* PUP POS E * 

* CORRECT MULTI-COMPONENT STRAIN GAGE BALANCE * 

* RECORD INGS FOR 1ST AND 2ND ORDER INTERACTIONS * 

* * 

* ASSUMPTION * 

* THE BALANCE RECORDINGS HAVE BEEN CONVERTED * 

* TO ENGINEERING UNITS. THAT IS, TUNNEL PRIME * 

* SENSITIVITIES HAVE ALREADY BEEN APPLIED * 

* * 

* LANGUAGE * 

* FORTRAN 2 OR A * 

* * 

* USAGE •* 

* DEFINE INPUT COMMON PARAMETERS * 

* CALL C INTRI FU»FZ»EZ,LIST,F» IERI * 

* * 

* DESCRIPTION OF INPUT CALLING SEQUENCE PARAMETERS * 

* FU UNCORRECTED COMPONENTS, ENGINEERING UNITS * 

* FZ CORRECT INITIAL LOADS, DETERMINED ITERATIVELY * 

* EZ 2ND ORDER INTERACTION DUE TO CORRECT INITIAL LOADS * 

* LIST PRINT OPTION TO DISPLAY THE PATTERN OF CONVERGENCE * 

* LIST=0 DO NOT PRINT COMPONENTS PER ITERATION * 

* OTHERWISE, LIST IS THE LOGICAL UNIT NUMBER * 

* * 

* DESCRIPTION OF OUTPUT CALLING SEQUENCE PARAMETERS * 

* F CCMPGNENTS CORRECTED FOR INTERACTIONS * 

* IER ERROR INDICATOR FOR INTERACTION CONVERGENCE * 

* I ER=0 INTERACTIONS CONVERGED * 

* IER= l INTERACTIONS OID NOT CONVERGE * 

* * 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

★ 

* 

* 

* 

* 

* 

* 

* 


DESCRIPTION OF INPUT COMMON PARAMETERS * 

IBAL BALANCE NAME IN DISPLAY CODE * 

EACH BALANCE HAS BEEN ASSIGNED A UNIQUE NAME * 

IDATE CALIBRATION DATE IN DISPLAY CODE * 

MCN T H/DAY/YEAR XX/YY/ZZ * 

KDATE CALIBRATION DATE EXPRESSED AS AN INTEGER * 

YEAR *10000+ MONTH* 100+0 AY ZZXXYY * 

M NUMBER OF BALANCE COMPONENTS PHYSICALLY DEFINED * 

M IS GREATER THAN 0, BUT LESS THAN OR EQUAL TO 6 * 

NAMEC ARRAY OF M COMPONENT NAMES IN A2 DISPLAY CODE * 

ALL MATRICES MUST BE ARRANGED ACCORDING TO NAMEC * 


ITASK INTEGER CODE SPECIFYING A TASK OR N1N2 TYPE BALANCE * 
ALL MATRICES MUST BE CONSISTENT WITH THE CALIBRATION* 
IORCR ORDER OF THE BALANCE CALIBRATION * 

ICRDR=0 NO INTERACTIONS * 

ICRDR=L 1ST ORDER INTERACTIONS ONLY * 

I CRDP =2 1ST AND 2ND ORDER INTERACTIONS * 

ITRNL OPTION TO TRANSLATE INTERACTIONS FOR INITIAL LOADS * 
I TRNL=0 DO NOT TRANSLATE FOR INITIAL LOADS * 

I TRNL= 1 DO TRANSLATE FOR INITIAL LOADS * 

IPLLS OPTION FOR ONE 2ND ORDER DISCONTINUOUS INTERACTION * 
I PLU S = 0 NO DISCONTINUOUS INTERACTION TERM * 

OTHERWISE, NAMEC! IPLUSI IS THE ACTING COMPONENT * 

MINUS INDEX IN C1IC2 TO ACCOMMODATE DISCONTINUITY * 

NCTE.ONE 2ND ORDER TERM CHANGES ONE COLUMN OF C1IC2 * 
NT P Y MAXIMUM NUMBER OF ITERATIONS ALLOWEO FOR CONVERGENCE* 
AN ERROR FLAG IS SET, IF NTRY IS INSUFFICIENT * 
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C 

C 

c 
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c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 
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c 

c 

c 

c 

c 
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♦ 

* 

* 

♦ 

* 

* 

* 

♦ 

* 

* 

* 

♦ 

* 

* 

* 

♦ 

* 

* 

♦ 

♦ 

* 

♦ 

♦ 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

+ 

* 

* 

* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

♦ 

* 

* 

* 

♦ 

* 


Cl I 

Cl IC 2 

CPOS 
CNEG 
PR CM 
CSENS 
ACCUR 


INVERSE OF NORMALIZED 1ST OROER INTERACTIONS WITH * 
MAIN DIAGONAL ELEMENTS OF 1. CONTAINS M* M ELEMENTS * 
PRODUCT OF C1I AND NORMALIZED 2ND ORDER INTERACTIONS* 


CONTAINS M*N ELEMENTS, WHERE N=M(M*L)/2 
ARRAY OF M POSITIVE CALIBRATION CONSTANTS 
ARRAY OF M NEGATIVE CALIBRATION CONSTANTS 
PERCENT ACCURACY REQUIRED FOR CONVERGENCE 
ARRAY OF M CALIBRATION PRIME SENSITIVITIES 
ARRAY OF M COMPONENTS REPRESENTING THE ACCURACY 
CF THE RECORDING SYSTEM, USED TO ESTABLISH THE 
INTERACTION CONVERGENCE CRITERIA. IT IS ASSUMED 
THAT ALL ELEMENTS OF ACCUR ARE GREATER THAN 0. 


BALANCE INTERACTION HISTORY FILE 

THE INPUT COMMON PARAMETERS RESIDE ON A BALANCE 
INTERACTION HISTORY FILE. THE FILE CONSISTS OF A 
PAIR OF RECORDS FOR EACH BALANCE. THE FIRST RECORD 
OF EACH PAIR CONTAINS THE ORIGINAL CALIBRATION MATRIX 
C t 16 2). THIS SUBROUTINE USES THE SECONO RECORD, 

WHICH CONTAINS THE INVERSELY DERIVED MATRICES C 1 I ( 36 ) 

AND C1IC2C126) 

REMARKS 

THIS SUBROUTINE IS DESIGNED IN SUCH A WAY THAT ALL 
COMPONENT AND CALIBRATION MATRICES COULD BE COLLAPSED 
TO CNLY THOSE COMPONENTS THAT ARE ACTUALLY HOOKED UP 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
GMECF EQUATES TWO MATRICES 
GMACF ADDS TWO MATRICES 
GM c BF SUBTRACTS TWO MATRICES 
GMPXF MULTIPLIES TWO MATRICES 
TESTF COMPARES TWO MATRICES 
GMSTF SETS A MATRIX EQUAL TO A SCALAR 

MATAS MULTIPLIES A MATRIX BY ITS TRANSPOSE AND STORES 

THE UPPER TRIANGLE OF THE PRODUCT 1-DIMENS I ON ALLY 

METHOD 

DETERMINING CORRECT COMPONENTS F IS VIEWED AS AN 
ITERATIVE SOLUTION TO THE FOLLOWING MATRIX EQUATION 

F - C1I X FU - C1IC2 X F2 
( M , 1 ) ( M , M) ( M, 1 ) <M,N) ( N , 1 ) 

WHERE F2 IS ALL PRODUCT COMBINATIONS OF F AND N*M(M+l)/2 


♦ 

* 

♦ 

* 

* 

♦ 

* 

* 

* 

♦ 

* 

♦ 

♦ 

* 

♦ 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 


LETTING EPSI 


C1IC2 X F2 


THE ITERATIVE TECHNIQUE FOLLOWS * 

♦ 


APPROXIMATION 


CL I 
Cl! 
Cl I 


FU 

FU - EPSI 
FU - EPS2 


C1I 
C1I 
Cl I 


RESULTS 

C FU - EPSI 
( FU - EPS2 
C FU - EPS3 


ERROR 

EPSI 

EPS2-EPSL 
EP S3-EPS2 


UNTIL A8S0LUTE<ePS<n-EPS< 1-1 ) ) < ABSOLUT E < ACCUR ) 
FOn ALL COMPONENTS 


* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 


******************************************************************* 


DIMENSION FU(6),FZ<6)»EZ(6) * F ( 6 ) 
DATA IZ/O/ 
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C 

c 

c 

c 
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c 
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c 

c 

c 

c 

c 
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c 

c 

c 
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INPUTS FPCM THE BALANCE INTERACTION HISTORY FILE 


COMMON/BAL/ 

IIBAL, ICATEtKCATE, M,NAMEC(6) , IT ASK, IOROR , I TRNL ♦ I P LUS , M I NUS, NTRY, 
2C1 II 36) ,C 1 IC2( 126 ) * CPOS <6) v CNEG(6 ) « PR CNT f CSENS ( 6) , ACC UR ( 6) 


WORKING STORAGE AREA AVAILABLE TO ALL SUBPROGRAMS 


COMMON/WORK/ 

LN,F2(21) * EPS! (6) * EP SO ( 6 ) , DELT A ( 6 ) » It J» ICNVG 


NO INTERACTIONS^ NOT NECESSARY TO CORRECT COMPONENTS 


I ER* 0 

IF(IORDR) 20,10,20 


•SET CCRRECTEO EQUAL TO UNCORRECTED AND RETURN 


10 CALL CMEQF(FU,F,M) 

GO TC 200 


CORRECT CCRFCNENTS FOR 1ST ORDER INTERACTIONS EXPLICITLY 


20 CALL GMPXFICl I,FU,F*M,M, 1) 


•THIS IS THE FIRST APPROXIMATION FOR 2ND ORDER CORRECTIONS 


I F ( ICRDR-1) 30,200,30 
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OPTION to translate axes, compensating for initial loaos 


30 miTRNLI 50,40,50 


TRANSLATION not necessary, initialize epsilon to zero 


40 CALL CNSTF!EPSO,0.,M) 

GO TC 60 


TRANSLATION NECESSARY, ADD INITIAL LOADS TO 1ST APPROXIMATION 


50 CALL GMACF(F,FZ,F,M) 


INITIALIZE EPSILON TO 2ND ORDER INTERACTION ON INITIAL LOADS 


CALL CMEQFI EZ,EPS0,M) 


OPTION TO HANDLE ONE DISCONTINUOUS 2ND ORDER INTERACTION TERM 


60 I F ( I PLUS ) 7C, ICO, 7C 


DETERMINE WHETHER TO USE POSITIVE OR NEGATIVE CALIBRATION 


70 I FI F I IPLLS) I 80,90,90 


SETUP TO USE INTERACTION TERM FROM NEGATIVE CALIBARTI ON 


29 



nnooooo o ooooooo onooaocio oooooooo onoonooo 


APPENDIX C - Continued 


80 CALL GMEQFICNEG,C1IC2(MINUS)»M) 

GC TC 100 
C 

c 

c 

c 

c 

c 

c 

C .SETUP TO USE INTERACTION TERM FROM POSITIVE CALIBRATION 

C 

C 

c 

90 CALL GMEQF(CP0S»C 1IC2I MINUS I , MI 


PRINT OPT I CN TO DISPLAY THE INTERACTION CONVERGENCE PATTERN 


100 IF(LIST) 11C, 130, 110 


ESTABLISH HE AO INGS OF ITERATION AND COMPONENT NAMES 


110 WRITEIL 1ST, 111) INAMECI J) , J=l,MI 

111 FORMAT ( 10H0 ITERATION, 8X, A2 , 5 ( 17X , A2 1 I 


PRINT FIRST APPROXIMATION OEEMEO ITERATION NUM8ER 0 


WRITE(LIST,120! I l, ( F C J I , J= 1 ,M I 
120 FORMAT! 16,6619.61 


CORRECT CCMPCNENTS FOR 2ND OROER INTERACTIONS ITERATIVELY 


130 DO 170 I = 1 ,MPY 


COMPUTE ALL PRODUCT COMBINATIONS OF THIS APPROXIMATION 
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APPENDIX C - Continued 


c 


CALL MATAS! F,M,F2,N) 


COMPITF 2ND ORDER INTERACTION DUE TO THIS APPROXIMATION 


CALL CMPXF(ClIC2,F2,EPSl,M,N,n 


COMPETE ERROR IN THIS APPROXIMATION OF CORRECTED COMPONENTS 


CALL CMSBF!EPSI,EPSO, DELTA, MI 


COMPETE NEXT APPROXIMATION OF THE CORRECTED COMPONENTS 


CALL GMSBFIF, DELTA, F,M) 


PRINT OPTION TO DISPLAY TOTAL LOADS PER ITERATION 


IF(LIST) 140,150,140 
140 WRITE(LIST,120) I , < F ( J I , J = 1 , M I 


DEMAND SIMULTANEOUS CONVERGENCE OF ALL COMPONENTS 


150 CALL TE STF I OELTA , ACCUR , K, ICNVGI 


DID INTERACTIONS CONVERGE TO WITHIN A PRESCRIBED ACCURACY 


I F ( ICNVGI 180,160,180 


NC, SAVE RESULTS OF THIS ITERATION AND TRY AGAIN OR 
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160 CALL GMECFIEPSI ,EPS0,M» 

.IF MAXIMUM TRYS EXCEEDED, SET ERROR FLAG AND RETURN 

170 CCNTINUE 

I E R= 1 

.TRANSLATE AXES BACK TO ORIGINAL SYSTEM 
180 IF(ITRNL) ISC *200 * 190 

.SUBTRACT INITIAL LOADS FROM TOTAL LOADS 

190 CALL CMSBF IF,FZ,F,M) 

• • • • •••••• ••••§•*••• •••••••• 

.RETURN COMPONENTS CORRECTED FOR INTERACTIONS 

200 RETLRN 

END 
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S UBR OUT IKE CTRNL(FZ,£Z) 

***4 c ******* 44 4* 4** **************************** **********4c********** 


* * 

* SUBROUTINE * 

* CTRM * 

* * 

* PURPOSE * 

* COMPUTE 2N0 CROER INTERACTION DUE TO INITIAL LOADS. * 

* PROVIDE INPUT TO SUBROUTINE CINTR FOR AXES TRANSLATION * 

* * 

* LANGUAGE * 

* FORTRAN 2 OR A * 

* * 

* USAGE * 

* DEFINE INPUT COMMON PARAMETERS * 

* CALL CTFNL(FZ,EZ) * 

* * 

* DESCRIPTION OF INPUT CALLING SEQUENCE PARAMETER * 

* FZ CORRECT INITIAL LOADS, DETERMINED ITERATIVELY * 

* * 

* DESCRIPTICN OF OUTPUT CALLING SEQUENCE PARAMETER * 

* EZ 2ND ORDER INTERACTION OUE TO CORRECT INITIAL LOAOS * 

* * 

* REMARKS * 

* THE INPUT CCMMON PARAMETERS RESIDE ON A BALANCE INTERACTION* 

* HISTORY FILE. SEE SUBROUTINE CINTR FOR A DESCRIPTION OF * 

* THESE PARAMETERS * 

* * 

* SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED * 

* GMECF EQUATES TWO MATRICES * 

* GMPXF MULTIPLIES TWO MATRICES * 

* GMSTF SETS A MATRIX EQUAL TO A SCALAR * 

* MATAS MULTIPLIES A MATRIX BY ITS TRANSPOSE AND STORES * 

* THE UPPER TRIANGLE IN I- D I MENS I ONAL SYMMETRIC FORM * 

* * 

* METHOD * 

* * 

* FZ = CLI X FUL - CIIC2 X F 2Z * 

* WHERE UNCORRECTED INITIAL LOAOS FU Z ARE NEVER * 

* REALLY KNOWN ANO F2Z IS ALL PRODUCT COMBINATIONS * 

* OF FZ. FOR THE PURPOSE OF TRANSLATING AXES IN * 

* SIBROUTINE CINTR IT IS SUFFICIENT TO KNOW FZ ANO ♦ 

* EZ = C1IC2 X F2Z BY DEFINITION * 

* * 


**********44**44******* ************************** ****************** 


DIMENSION FZ(6),EZ<6) 


.INPUTS FRCM the BALANCE INTERACTION HISTORY FILE 


COMMON / 8 At / 

II BAL , I DATE ,KC ATE, M ,NAMEC (6) , I T ASK 1 1 OR DR , I T RNL , I PLUS , M I NUS , NTRY , 
2C1I (36 ) , Cl IC2( 126 ),CP0S<6J , CNEG ( 6 ) , PR CNT, C SEN S ( 6 ) ,ACCUR(6) 

C 

C 
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APPENDIX C - Continued 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


WORKING STORAGE AREA AVAILABLE TO ALL SUBPROGRAMS 


COMMON/WOPK/N *F2Z(2l) 


TRANSLATION ONLY NECESSARY FOR 2NO ORDER INTERACTIONS 


I F ( I OP DR— 2 ) 10,20, 10 


OTHERWISE, SET INITIAL EPSILON TO ZERO 


10 CALL GMSTF(EZ,0. f M) 

GO TC ICO 


OPTION TO HANOLE ONE DISCONTINUOUS 2ND ORDER INTERACTION TERM 


20 IF(IPLUS) 30,60,30 


DETERMINE WHETHER TO USE POSITIVE OR NEGATIVE CALIBRATION ! 


30 IF(FZUPLIS) ) 40,50,50 


SETUP TO USE INTERACTION TERM FROM NEGATIVE CALIBRATION^ 

* 


40 CALL GM EOF t C NEC t Cl IC2J MINUS I »M1 

GC TC 60 


SETUP TO USE INTERACTION TERM FROM POSITIVE CALIBRATION^ 


50 CALL GMEQF(CP0S,C1IC2(MINUS),M> 
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APPENDIX C - Continued 


COMPUTE ALL PROOUCT COMBINATIONS OF CORRECT INITIAL LOADS 


60 CALL MATAS (FZ,M,F2Z,N) 


C 

C 


C 

c 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


COMPUTE 2NC CRDER INTERACTION OUE TO CORRECT INITIAL LOADS 


CALL GMPXF (C1IC2,F2Z, EZ, M,N, 1) 

100 RETURN 
END 

SUBROUTINE GMEQF { A , R , MN ) 

********** ****************************** **44 t * **************** 


* * 

♦ SUBROUTINE * 

♦ GMEQF * 

♦ * 

* PURPOSE * 

* EQUATE CNE GENERAL MATRIX TO ANOTHER GENERAL MATRIX * 

* * 

* LANGUAGE * 

* FORTRAN 2 OR A * 

* * 

* USAGE * 

* CALL GME CF ( A * R t MN ) * 

* * 

* DESCRIPTION OF PARAMETERS * 

* A INPUT MATRIX NAME * 

* R OUTPUT MATRIX NAME * 

* MN INPUT NUMBER OF ELEMENTS IN MATRIX A OR R * 

* * 

♦ REMARKS * 

* 1. THE ELEMENTS OF MATRIX A ARE NOT CHANGED. * 


* 2. THE USER IS CAUTIONED, IF MATRICES A AND R ARE NOT FLOATING * 

* POINT. FOR EXAMPLE, TWO INTEGER TO ONE FLOATING POINT WORD.* 

* 3. SUBROUTINE GMEQF CAN BE USED TO MANIPULATE MATRIX COLUMNS. * 


* FOP EXAMPLE, SET MATRIX R<M,i) EQUAL TO THE JTH COLUMN OF * 

* MATRIX A ( M , N ) BY CALL GMEQF (A( 1 , J ) , R , M ) . * 

* 4. SUBROUTINE GMEQF CANNOT BE EASILY USED TO MANIPULATE ROWS * 

* DUE TO TEE FACT THAT THE ELEMENTS OF A ROW ARE NOT HELD * 

* CONS ECLTIVELY IN CORE STORAGE. * 

* * 

* SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIREO * 

* NONE * 

* * 

* METHOC * 

* EACH ELEMENT OF MATRIX R IS SET EQUAL TO * 

* THE CCFRESPCNCING ELEMENT OF MATRIX A * 

* * 

* RUJ)*A(IJ) FOR IJ=I,2,...,MN * 

* * 


******************************************************************* 


35 



o o o 


APPENDIX C - Continued 


EQUATE MATRICES 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DIMENSION A ( l ) ,R ( 1 ) 

DO 10 IJMtMN 
R ( l J ) * A( I J ) 

10 CONTINUE 
RETURN 
END 

SUBROUTINE GRACE! A , B# Rt MN ) 

************** ***************************************************** 


* + 

* SUBROUTINE * 

★ GMACf * 

♦ * 

* PURPOSE * 

* ADC TWC GENERAL MATRICES TO FORM RESULTANT GENERAL MATRIX * 

* * 

* LANGUAGE * 

* FORTRAN 2 OR A * 

* * 

* USAGE * 

* CALL C MA CF ( A » B » R * MN ) * 

* * 

* DESCRIPTICN OF PARAMETERS * 

* A INPUT FIRST MATRIX NAME * 

* B INPUT SECOND MATRIX NAME * 

* R OUTPUT MATRIX NAME ♦ 

* MN INPUT NUMBER OF ELEMENTS IN MATRIX A,B,OR R * 

* * 

* REMARKS * 

* MATRICES A t B * AND R MUST BE FLOATING POINT * 

* MATRICES A«e« AND R MAY BE THE SAME LOCATIONS * 

* OTHERWISE, THE ELEMENTS OF MATRICES A f B ARE NOT CHANGED * 

* * 

* SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED * 

* NONE * 

* * 

* METHOD * 

* EACH ELEMENT OF MATRIX A IS ADDED TO THE CORRESPONDING * 

* ELEMENT OF MATRIX B AND THE RESULT IS PLACED IN THE * 

* CORRESPONDING ELEMENT OF MATRIX R * 

* * 

♦ R ( I J ) =A( IJ)+B( I J ) FOR IJ*l f 2i...,MN * 

♦ * 


******************************************************************* 

ADD MATRICES 

DIMENSION A<l)fB(l),R(l) 

DO 10 IJMtMN 
RUJ) = A( IJI + e( IJ) 

10 CONTINUE 
RETURN 
END 


36 



APPENDIX C - Continued 


SUBROUTINE GMSBF ( A , B , R , MN) 

C 

c ************************* ***************************** ************* 


c 

♦ 



♦ 

c 

♦ 

SUBROUTINE 


* 

c 

* 

GMSeF 


* 

c 

* 



* 

c 

* 

PURFOSE 


* 

c 

* 

SUBTRACT ONE GENERAL MATRIX FROM ANOTHER 


♦ 

c 

* 

TO FORM t RESULTANT GENERAL MATRIX 


♦ 

c 

* 



♦ 

c 

♦ 

LANGUAGE 


* 

c 

* 

FORTRAN 2 OR 4 


* 

c 

♦ 



* 

c 

* 

USAGE 


* 

c 

* 

CALL GMS 6F ( A f B» R * MN ) 


♦ 

c 

* 



* 

c 

* 

DESCR I FT I CN OF PARAMETERS 


* 

c 

* 

A INPUT NAME OF FIRST MATRIX 


* 

c 

♦ 

B INPUT NAME OF SECOND MATRIX 


* 

c 

* 

R OUTPUT MATRIX NAME 


♦ 

c 

* 

MN INPUT NUMBER OF ELEMENTS IN MATRIX A,B,OR 

R 

* 

c 

* 



* 

c 

* 

REMARKS 


* 

c 

* 

MATRICES A » B » AND R MUST BE FLOATING POINT 


* 

c 

♦ 

MATRICES A,8» AND R MAY BE THE SAME LOCATIONS 


♦ 

c 

♦ 

OTHERWISE, TEE ELEMENTS OF MATRICES A,B ARE NOT 

CHANGED 

* 

c 

♦ 



* 

c 

♦ 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 


* 

c 

♦ 

NONE 


* 

c 

* 



* 

c 

* 

METHOD 


* 

c 

♦ 

EACH ELEMENT OF MATRIX B IS SUBTRACTED FROM THE 


* 

c 

♦ 

COPPESPCNCING ELEMENT OF MATRIX A ANO THE RESULT 


* 

c 

* 

IS PLACED IN THE CORRESPONDING ELEMENT OF MATRIX 

ft 

* 

c 

♦ 



* 

c 

* 

R(IJ)=A( I JJ-BI IJI FOR I J=l,2, ...,MN 


* 

c 

* 



* 


c ******************************************************************* 

c 

C SUBTRACT MATRICES 

C 

0 IMENS ION All), Bill, R(l) 

00 10 IJ=1,MN 

r ( i j ) = a ( u)-e< m 
10 CONTINUE 
RETURN 
END 
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APPENDIX C - Continued 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE GMPXFt A»B,R»MfN,L) 


***************** ***************************** ********************* 

♦ 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


SUBROUTINE 
GMP X F 

PURPOSE 

MULTIPLY TWO GENERAL MATRICES 
TO FORM A RESULTANT GENERAL MATRIX 


LANGUAGE 

FORTRAN 


2 OR 4 


USAGE 

CALL GMPXF { A, B,R ,M f N, L) 
DESCR IPTICN OF PARAMETERS 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 

METHOD * 

TEE M PY N MATRIX A IS POSTMULT I PLI ED BY THE N BY L * 

MATRIX 0 AND THE RESULT IS STORED IN THE M BY L MATRIX R. * 

* 

* 
* 
* 


A 

B 

R 

M 

N 

L 


INPUT FIRST MATRIX NAME 
INPUT SECOND MATRIX NAME 
OUTPUT MATRIX NAME 

INPUT NUMBER OF ROWS IN MATRIX A OR R 
INPUT NUMBER OF COLUMNS IN A AND ROWS IN B 
INPLT NUMBER OF COLUMNS IN MATRIX 0 OR R 


REMARKS 

ALL MATRICES MUST BE STORED IN FLOATING POINT 

A ANC B MUST BE CONFORMABLE FOR MATRIX MULTIPLICATION 

A ANC B MAY BE THE SAME MATRIX IF IT IS SQUARE 

MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRIX A OR B 

THE ELEMENTS OF MATRICES A AND B ARE NOT CHANGED 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NCNE 


FCR A GIVEN ROW I AND COLUMN J, 
R<ItJ)-THE SUMMATION FROM K=lt2'.< 
CF THE PRODUCTS MI|K)*B(K»J) 


****************************************************** ************* 


MULTIPLY MATRICES 

DIMENSION Ml)fB(i) *R( l) 
I R=0 
IK®— N 

DO 30 K=1 ,L 
I K= I K + N 
DO 20 J® 1 ? M 
I R= I R ♦ l 
J I 3 J— M 
IB® IK 
R(IR)=0. 

DO 10 1= 1 1 N 
J 1 = J I * M 
I B® I B+ 1 

RUR>*R<m + MJM*Bt IB) 
10 CONTINUE 
20 CONTINUE 
30 CONTINUE 
RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE T ESTF ( A , B, MN, LE ) 

******************************************************************* 


* * 

* SUBROUTINE * 

* TESTE * 

* * 

* PURPOSE * 

* TEST Tl-E ABSOLUTE VALUE OF EACH ELEMENT OF MATRIX A TO * 

* DETERMINE IF IT IS LESS THAN OR EQUAL TO THE CORRESPONDING * 

* ELEMENT CF MATRIX B * 

* * 

* LANGUAGE * 

* FORTRAN 2 OR 4 * 

* * 

♦ USAGE * 

♦ CALL TESTFI A,E,MN,LE) * 

♦ * 

* DESCRIPTICN OF PARAMETERS * 

* A INPUT FIRST MATRIX NAME * 

* B INPUT SECOND MATRIX NAME * 

* MN INPUT NUMBER OF ELEMENTS IN MATRIX A OR B * 

* LE OUTPUT COMPARISON OF MATRICES A AND B * 

* * 

* REMARKS * 

* LE=C IF THE ABSOLUTE VALUE OF AT LEAST ONE ELEMENT IN * 

* MATRIX A IS GREATER THAN THE VALUE OF THE CORRESPONDING * 

* ELEMENT IN MATRIX B. OTHERWISE, LE=I * 

* * 

* SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED * 

* NONE * 

* * 

* METHOD * 

* IF I A ( I J ) 1 LESS THAN OR EQUAL TO B(IJ) * 

* FOR ALL I J=l,2, ... , MN THEN LE=1 * 

* OTHERWISE, LE*Q * 

* * 


********** 44 ******************************************************* 

COMPARE MATRICES 
DIMENSION A ( 1 ) f 8 ( I ) 


L E= 0 

DO 10 I J = 1 » MN 

I F( ABS( A ( IJ) 1— B I IJ) ) 10,10,20 
10 CONTINUE 
L E= 1 

20 RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE C m ST F < R , S f MN ) 

*********44 4 ***************** ********** * ***** ********************** 


♦ * 

* SUBROUTINE * 

* GNSTF * 

* * 

* PURPOSE * 

* SET ALL ELEMENTS OF A GENERAL MATRIX EQUAL TO A SCALAR * 

* * 

* LANGUAGE * 

* FORTRAN 2 OR A * 

* * 

* USAGE * 

* CALL GM$TF(R»Sf MN) * 

* * 

* DESCRIPTICN OF PARAMETERS * 

* R OUTPUT MATRIX NAME * 

* S INPUT SCALAR CONSTANT * 

* MN INPUT NUMBER OF ELEMENTS IN MATRIX R * 

* * 

* REMARKS * 

* ALL VARIABLES SHOULD BE FLOATING POINT * 

* * 

* SUBROUTINES ANO FUNCTION SUBPROGRAMS REQUIRED * 

* NONE * 

* * 

* METHOD * 

* SET EACH ELEMENT OF MATRIX R EQUAL TO THE SCALAR S * 

* ♦ 

* R < I J ) “S FOR ALL I J* 1 1 2 • • * . * MN * 

* * 


***************** ****************************************** ******** 

SET EACH ELEMENT OF MATRIX R EQUAL TO THE SCALAR S 

DIMENSION F C 1 1 
DO 10 I J = 1 1 MN 
R ( I Jl *S 
10 CONTINUE 
RETURN 
END 
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APPENDIX C - Continued 


SUBROUTINE MATAS l A,M,R,NI 

******************************************************************* 
* * 

* PURPOSE * 

* POSTMU IT I PLY A COLUMN MATRIX BY ITS TRANSPOSE AND STORE THE* 

* UPPER TRIANGLE OF THE RESULTANT MATRIX IN SYMMETRIC FORM * 

* * 

* LANGUAGE * 

* FORTRAN 2 OR A * 

* * 

* USAGE * 

* CALL MATASI A,M,R,N) * 

* * 

* DESCRIPTION OF PARAMETERS * 

* A INPUT MATRIX NAME * 

* M INPUT NUMBER OF ELEMENTS IN MATRIX A * 

* R OUTPUT MATRIX NAME * 

* N OUTPUT NUMBER OF ELEMENTS IN MATRIX R * 


c 

♦ 

EXAMPLE l-l 1- -1 

1- 





-1 

1--L 

* 

c 

♦ 

1M X IN A P R Y SI = 

1NN 

NA 

NP 

NR 

NY 

NS 1 

= INN 1 * R 

* 

c 

* 

1 A 1 1- -l 

1 

AA 

AP 

AR 

AY 

AS 1 

1 NA 1 

♦ 

c 

* 

1P1 ( L , 6 ) 



PP 

PR 

PY 

PS1 

1NP 1 

* 

c 

* 

1*1 




RR 

RY 

RSI 

INR 1 

* 

c 

* 

1Y1 





YY 

YS 1 

1NY 1 

* 

c 

* 

1S1 






SSI 

1 NS 1 

* 

c 

* 

1-1 1-1 






- 1 

1 AA 1 

♦ 

c 

* 

(6,1) INI 





(6 

*6) 

1 AP l 

* 

c 

* 

1A1 







1 AR I 

* 

c 

♦ 

WHERE A = 1PI 







1 AY 1 

* 

c 

♦ 

1R1 







IASI 

* 

c 

* 

1Y1 







1PP1 

* 

c 

* 

IS L 







1 PR 1 

♦ 

c 

* 

1-1 







1PY1 

♦ 

c 

* 








1PS1 

* 

c 

♦ 

M=6 , N=2 1 







1RR1 

* 

c 

* 








1RY 1 

* 

c 

♦ 








1RSI 

* 

c 

* 








L YY 1 

* 

c 

* 








1YSI 

* 

c 

* 








1SS1 

* 

c 

♦ 








1—1 

* 

c 

♦ 








(21,1) 

* 

c 

* 









* 

c 

♦ 

REMARK? 








* 


l- THE RESULTANT NUMBER OF ELEMENTS IN MATRIX R IS N=M(M*ll/2 * 

2. FOR COMPUTER EFFICIENCY, MATRIX A IS RESTRICTED TO 1 COLUMN* 

3. MATRICES A AND R CANNOT SHARE THE SAME LOCATIONS * 

A. THE ELEMENTS OF MATRIX A ARE NOT CHANGEO * 

5. MATRIX R REPRESENTS ALL PRODUCT COMBINATIONS OF M ELEMENTS * 

* 

FUNCTIONS AND SUBPROGRAMS REQUIRED * 

NONE * 


METHOO * 

ANY MATRIX AIM, LI TIMES ITS TRANSPOSE ATI L , MI RESULTS IN A * 
SYMMETRIC MATRIX R I M , M ) . THIS SUBROUTINE HAS RESTRICTED * 
L T C 1. THE UPPER AND LOWER TRIANGLES OF MATRIX R ARE * 

IMAGES OF ONE ANOTHER. CERTAIN APPLICATIONS REQUIRE * 

ONLY THE UPPER OR LOWER TRIANGLE STORED l-DI MEN SI ONALLY . * 

THE ABOVE EXAMPLE DEMONSTRATES THE 1-DIMENSIONAL ORDERING * 
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APPENDIX C - Concluded 


PERFORM THE MATRIX OPERATION 

DIMENSION M 1 ) » R ( 1 ) 

N=0 

DO 20 1 = 1, M 
DO 10 J=I,P 
N=N+1 

R (N ) - A ( I )*AU) 

10 CONTINUE 
20 CONTINUE 
RETURN 
END 
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